Extrinsic noise driven phenotype switching in a self-regulating gene 
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Due to inherent noise in intracellular networks cellular decisions can be random, so genetically 
identical cells can display difTerent phenotypic behavior even in identical environments. Most previ- 
ous work in understanding the decision-making process has focused on the role of intrinsic noise in 
these systems. Yet, especially in the high copy-number regime, extrinsic noise has been shown to be 
much more significant. Here, using a prototypical example of a bistable self-regulating gene model, 
we develop a theoretical framework describing the combined effect of intrinsic and extrinsic noise on 
the dynamics of stochastic genetic switches. Employing our theory and Monte Carlo simulations, 
we show that extrinsic noise not only significantly alters the lifetimes of the phenotypic states, but 
can induce bistability in unexpected regions of parameter space, and may fundamentally change 
the escape mechanism. These results have implications for interpreting experimentally observed 
heterogeneity in cellular populations and for stochastic modeling of cellular decision processes. 
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Introduction. Noise-driven switching between coexist- 
ing metastable states plays a key role in many systems in 
physics, chemistry, and biology 1-4 . Besides thermal or 
intrinsic noise (IN) that drives switching [5 , such systems 
often experience extrinsic or environmental noise (EN) 
from the noisy environment or from being coupled to an- 
other fluctuating system [Ij. Noise-driven escape from 
a metastable state while under the influence of EN has 
been previously studied in the context of population bi- 
ology and population genetics (see e.g. [5HH]). Here, e.g., 
it has been shown that delta-correlated as well as col- 
ored EN can drastically decrease the population's mean 
extinction time [71IH]. Moreover, recently there has been 
a large effort to predict the onset of EN-driven critical 
transitions and regime shifts in ecosystems, see e.g. [1]. 

In cellular biology, most studies of gene expression 
dynamics, including our own treatments |10j . have fo- 
cused on the role of IN (reviewed in [TI]). Recently, 
however, gene expression under EN has also come un- 
der study [12l [13], where EN has been experimentally 
conflrmed to be one of the dominant sources of variation 
in protein copy number, particularly above copy numbers 
of 0(10) [l3]. In studies of genetic switches, EN has been 
shown to induce bistability ^14. ,15) , vary the distribution 
tails pS] and modify switching times [16]. Yet, previous 
studies have not provided fundamental insight as to the 
interplay between IN and EN in the switching process, 
i.e., how the mean switching times (MSTs) and switching 
paths deviate according to EN strength, correlation time 
and statistics. Elucidating the relationship between IN 
and EN during switching is crucial to understanding how 
EN affects population heterogeneity in bistable systems, 
which is of importance when studying, e.g., bet-hedging 
strategies like bacterial persistence [T7]. 

In this Letter we study the contributions of IN and 
EN to noise-driven switching in a simple self-regulating 



genetic circuit with positive feedback. Employing a semi- 
classical theory we perform a systematic study of the ef- 
fect of EN statistics, magnitude and correlation time, on 
the switch's stochastic dynamics. In particular, we derive 
expressions for the MSTs as functions of the EN strength 
and correlation time, and also study how EN can induce 
bistability in an otherwise monostable system. All ana- 
lytical results are corroborated by extensive Monte-Carlo 
(MC) simulations. Our main conclusion is that EN cor- 
relation time plays a significant role in determining both 
the stability of the metastable state and the mechanism 
of escape. This strongly indicates that in biological sys- 
tems, where the correlation time is thought to be long, 
phenotype switching may be driven primarily by EN. 

Model. Our analysis relies on the model of a self- 
regulating gene (SRG), with positive feedback due to the 
production rate depending on the state. Let nit) be the 
protein copy number and N be the protein abundance 
in the hi state. Proteins are produced at a rate f{n), 
which is any Hill- like function, and decay with rate 1. 
The mean protein concentration x{t) — n{t)/N satisfies 

x^f{x)-x. (1) 

For simplicity we take f{x) ~ aQ+{l—aQ)0(x—XQ), where 
0{x) is the Heaviside step function, and ao < 2:0 < 1. 
Eq. ([!]) leads to a bistable system with three fixed points 
xi < X2 < X3, where xi = and 2:3 = 1 are attracting 
fixed points of the low and hi states respectively, while 
X2 = xq is repelling. Typically, ao <C 1 so X3 ^ xi. 

To account for IN, we employ the master equation for 
Pn(t) - the probability to find n proteins at time t: 

Pn = f{n - l)P„_i + (n + l)P„+i - [fin) + n]P^. (2) 

For simplicity we focus on the weak-noise regime 1— <C 
1, where (without loss of generality) the "switching bar- 
rier" between the hi and low states is small. In this 
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regime Eq. ^ is accurately approximated [TH| by the 
following Fokker-Planck equation (FPE) for the proba- 
bility P{x,t) to find concentration x at time t [3]: 

dtP = -dA[f{x)-x]P} + l/{2N)dl{[f{x) + x]P}. (3) 

Starting from the vicinity of the hi state, the sys- 
tem rapidly forms a quasi-stationary distribution (QSD) 
about the hi state, which slowly leaks through the un- 
stable point X — xq [T51 - BD] . In general, the metastable 
state decays as P{x,t) ~ 7r(x)e~*/'^ where t:{x) is the 
QSD and r is the MST. Employing the WKB ansatz 
7r(a;) - e-^^(=") for the QSD, where S{x) is called the 
action and p{x) = S'{x) is called the momentum |19j . 
Eq. ([3]) gives rise to a stationary Hamilton- Jacobi equa- 
tion (HJE) H{x,px) = with Hamiltonian 

H{x,p^) = p^[f{x) -x] + {pl/2)[f{x) + x]. (4) 

Switching occurs along the zero-energy trajectory 
p^{x) = -2[f{x) - x]/[fix) + x] of 0. For xq < x < 1, 
Px{x) = —2(1 — x)/{\ + x), which for 1 — <C 1 sat- 
isfies ^ 1. This yields S{x) = px{x')dx' = 
2\x — 21n(l -|- a;)], and the QSD around x ^ 1: 7r(x) ~ 
^-N[S{x)-s{i)] standard deviation CTi„ = N''^/'^. 
Therefore, since Thi^iow ^ 7r(xo)~^ [ISlUn], we have |5T] 

\nTh^^io^c^N[S{xoyS{l)]c^{N/2){l~x^f = ASo. (5) 

which is applicable as long as tTi„ = N~^/^ <C I — Xq. 

Next, we incorporate EN in the form of one or more 
fluctuating parameters. We assume that cell-to-cell vari- 
ability in transcription and translation rates causes the 
protein production rate to fluctuate. In the hi state the 
production rate then becomes ai{t) = 1 + £,{t), where 
^(t) is fluctuating with flnite correlation time. As we are 
interested in the hi — >■ low transition we ignore fluctua- 
tions in ao 1- We take ^(i) to be Ornstein-Uhlenbeck 
(OU) noise ^ : positively correlated Gaussian noise with 
zero mean, variance and correlation time Tc, satisfy- 
ing = crLe"'*"*''/^'^- The OU process satisfies 
the following Langevin equation 



(6) 



where rj is white Gaussian noise, {ri{t)rj{t')) = 5{t—t') [22] . 
Here, and Tc are characteristic of the environment 
and the cell's regulatory network and are generally un- 
known. Non-Gaussian statistics for EN have also been 
proposed |15j . but further theoretical and experimental 
work is needed to uncover the source and form of EN. 

To study the interplay between IN and EN, we com- 
bine Eq. ^ with the underlying IN dynamics [Eq. ([3|]. 
Defining the fluctuating production rate /(a;,^) — ao -I- 
(1 — ofQ + ~ ^0)7 drift term A{x, = /(x, ^) — x, 

diffusion coefficient B{x,^) = f{x,£,) + x, and the EN 
and IN variance ratio, V = (Jex/'^l 
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FIG. 1: (Color online) Thi^iow versus the relative strength for 
white Tc = 10"^ (a) and long-correlated Tc = 10^ (b-|-c) EN. 
MC simulations (the symbols) for noise in the birth (o) and 
death (x) terms (that are indistinguishable), are compared 
with theory (lines): Eq. ^ (a) and Eq. ([l2| (b-l-c). Here 
N = 5000, ao = 0.01, and xq = 0.93 (a-l-b) or xo = 0.915 (c). 



a 2-D FPE for the joint probability P{x, ^, t) to find con- 
centration X and noise magnitude 1^ at time t j231 124] : 



d^p^~d^{A{x,OP} 

+l/{2N)dl{B{x,0P} 



di mrc)P} 

l/i2N)dUi2V/r,)P}. (7) 



Employing the WKB ansatz P{x, £_) - e"^'^^"''^) for the 
QSD, Eq. yields a HJE: H{x,^,p^,p^) ^p^A{x,^)- 
£,P^/tc + {p^/2)B{x,^)+p'^V/tc = 0, with momenta = 
dxS and p^ = d^S. The HJE can be solved by considering 
the Hamilton equations Xi — dp-H and pi — —dxiH: 



x = A + pxB , px = 



-Px[dxA+{px/2)dxB] 



(8) 



where we have combined the equations for ^ and to a 
single equation for ^ and kept terms up to ©(p^) <C 1. 

Eqs. (H can be solved numerically for generic noise, 
which yields the corresponding action function 5(a;,^) = 
/ Px{x,£,)dx + p^{x, ^)d^, and QSD. Analytical progress 
can be made in two limits: short-correlated white noise 
Tc ^ 1, and long-correlated adiabatic noise Tc >• 1. 

For white EN, we neglect ^ in the third of Eqs. (jsj j8], 
which yields ^ ~ 2pxVTc. Substituting ^ into the first 
of Eqs. ([S]), we find for x > xq: x = f{x) ~ x + 
2pxVTc+Px[fix)+x + 2pxVTc], which originates from an 
effective white-noise Hamiltonian: H ~ Px[f{x) — x] + 
{Px/2)[f{x) -I- X -I- 2Vtc], where we have neglected 0{px) 
terms. Solving H = 0, we find Px{x) = —2(1 — x)/{l + 
x + 2Vtc), which yields the MST in the whitc-EN regime 



(9) 



Na^^, we obtain 



Eq. ([9]) is confirmed by MC simulations [25], see Figs. 
[T}|]2] In Fig. [2] and below, fi denotes the QSD's average. 

Now, to deal with long-correlated EN, we note that 
when Tc ^ 1, during the rare fluctuation that takes the 
system from the hi to the low state, the system sam- 
ples an almost constant value of the noise ^ = ■Co [H]- 
For a constant ^q, the hi flxed point becomes 1 -I- ^o- 
The optimal value of is found by minimizing the cost 
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FIG. 2: (Color online) Thi^iow as function of Tc for various EN 
strengths: (Jex/ 1^- = 0.01 (top left), CTex/z^ = 0.1 (top right), 
and a^x j^i = 0.2 (bottom left). The lines are the analytical 
predictions of Eq. Q (solid) and Eq. ( 12 \ (dotted). The lower 



right panel shows rj"^ versus the relative EN strength. Nu 
merical results (symbols) agree well with Eq. ( 13 1 (line) multi 
plied by 0.72. Inset shows the clear scaling of r°^* ~ (ct^i^/m)^ 



of switching given noise realization liiT/ji_j./oui(^o) — 
(Af/2)(1 — + ^o)^, against the (absolute value of the) 
statistical weight of ^, NQ^I(2V^. By doing so, we find 
^opt ^ _(^^ _ x^)v/ll + V), where < 1 - xq as 

expected. Plugging into Thi^iowi^o) we find [29] 



(10) 



For strong EN, > 1, Eq. ([lO| holds when ASq > V^, 
which can only be satisfied when aex ^ ^ ^ xq- 

What happens when aex ^ 1 — a^o 3> cTm? Here, IN 
can be neglected, and the MST turns out to be dominated 
solely by EN. Namely, the MST can be approximated by 
the mean first passage time T{x) it takes the OU process 
to reach position x starting from x = at i = 0. Using 
Eq. (|6|, T{x) is governed by the following equation [3]: 



{aiJre)T"{x)-{x/Te)T'{x)^^l, 



(11) 



with boundary conditions T{0) — and r'(oo) — 0, 
whose solution is T(x) = Tcip{x,aex)- Here, ip{x,(Tex) = 
(V2)Erfi(z) - z^aFa [{1, !},{§, 2}, , where z = 

x/iV2aex), Erfi(0) = 2/^nJ^ey"dy, and 2f2({}, {}, 2:) is 
the generahzed hypergeometric function. The MST is ob- 
tained by plugging x = I-xq: Thi^iow ^ T{l-xo) ~ 7rc, 
with 7 = 0{\) for a^x — 1 — Xq. 



This analysis gives rise to a correction in Eq. ( 10 ) 



for the MST in the adiabatic regime Tc ^ 1. Since, 



In Tc at aex ^ ^ — xq, and In t/i; 



at aex = 0, by defining A = In Tc, Eq. (10) becomes 



lnT;,,^,„^ ~ A + (A^o - A)(l + Vy 



(12) 



Eq. ( 12 1 compares well with MC numerics, see Figs. [ljl|2| 
As can be seen in Fig. [2] for given EN strength aex 
there exists an optimal EN correlation time Tc for which 




FIG. 3: (Color online) (a+b) The location of the stochastic 
fixed points of Eq. ([l]) with f{x) given by (14 1, for EN with 



10'' and aex/fJ. = 0.01 (left) and a^^/^ = 0.05 (right). 
The solid line shows the deterministic fixed points, (c-e) The 
steady-state 2-D PDFs of finding protein number n and in- 
stantaneous EN magnitude for a^^j [i, = 0.05 and xq — 0.42, 
0.48, and 0.56, respectively. Here N = 300, and ao = 0.05. 



the MST is minimal. In order to calculate r"^* we add 
the white- and adiabatic-noise contributions [Eqs. ^ and 
( 12 )] for the MST, and differentiate the result with re- 
spect to Te- For 1 — Xq ^ 1, wc find 



opt 



(13) 



whose dependence on aex is confirmed by Fig. [2] 
Noise in the degradation rate. We now consider the case 
where the degradation rate is fluctuating as l-l-^(t). Here, 
the corresponding FPE is given by Eq. ([t]) with A{x, ^) = 
fix) - xil + C) and B{x, ^ - f{x) + x{l + C). 

In the white EN regime the optimal path for switching 
at X > Xq becomes Px{x) = —2(1 — x)/{l + x + 2x'^Vtc)- 
This yields a MST that coincides with Eq. (|9|, see 
Fig. [1] since the x'^ factor in the denominator of Px{x) 
approximately equals 1 along the integration regime 
1 — Xq ^ 1. For adiabatic EN, the hi fixed point be- 
comes 2:3 = (1 + ^o)^^- This again yields after some 
algebra ~ [l ~ xq)V / [l + V), which coincides up to 
a minus sign with when the production rate is fiuc- 
tuating. Therefore, we recover Eq. (10 1, see Fig. [l] 
Noise-induced bistability. To study how EN affects bista- 
bility in the SRG model we worked with a modified 
production rate f{x) that allowed bidirectional transi- 
tions between the low and hi states with MSTs that 
were reachable using MC simulations. Instead of a step- 
function for f{x) in Eq. (jlj, we took 



f{x) = ao + (1 - ao)x'^/ (x^ + Xq), 



(14) 



where ^ 1- Given ao this system is bistable over a 
range of xq values. We are interested in how the range of 
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bistability varies under the influence of EN, and also how 
the latter affects the MSTs and the steady state proba- 
bility distribution functions (PDFs). To answer these 
questions we ran MC simulations with a degradation rate 
1 + C(^) in the adiabatic limit with Tc — 10'^. 

To determine the effect of the EN on the bistability 
range, we calculated the PDFs at various xq values from 
long-time simulations and extracted the position(s) of the 
sole maximum (monostable) or the two maxima sepa- 
rated by a minimum (bistable). These values we inter- 
preted as stochastic equivalents to the deterministic fixed 
points. For very weak EN, the stochastic and determin- 
istic fixed points generally agree with only small devia- 
tions. Yet, as the EN strength increases the locations of 
the stochastic fixed points undergo a dramatic departure 
from their deterministic locations. The example shown 
in the upper panels of Fig. [3] demonstrates that even for a 
modest EN strength of (Jex/lJ- = 0.05, the range of xq over 
which the system is bistable has greatly increased. This 
effect becomes more pronounced as the EN strength fur- 
ther increases. For Cex/M = 0.2 the system was bistable 
over the entire range of xq sampled (0.3-0.7). 

To further investigate the change in switching behav- 
ior, we calculated the 2-D PDFs of finding protein num- 
ber n and instantaneous fluctuation magnitude ^i. The 
lower panels of Fig. [3] show that, for strong EN, has a 
direct impact on the state of the system. Fig. [sj^d) shows 
a case where the system is deterministically bistable. 
Here, when is relatively weak the system undergoes 
noise-driven switching as expected. However, when the 
degradation reaction is sampling the highest rates, the 
system exists only in the low state and vice versa. When 
the EN drives the degradation rate to one of its extremes 
the system switches deterministically to the appropri- 
ate stable state. This effect appears in the 2-D PDFs 
as two alternate switching paths: when > 0, there 
is a hi ^ low pathway for leakage of probability, but 
when < there is a separate low — > hi leakage path. 
Thus, the system's bistability is not only a consequence 
of stochastic switching between states, but also of EN 
driving the system between different regions of parame- 
ter space with alternate fixed point configurations. 

Fig. [3](c-|-e) show the case where the system is deter- 
ministically monostable. When is low one can see that 
the system behaves as though it has a single fixed point. 
However, when a large fluctuation occurs in the correct 
direction it can shift the system into a region of param- 
eter space that is bistable; the fluctuations induce bista- 
bility in the system. This effect gives rise to the greatly 
increased bistability range observed in the simulations. 

Finally, we calculated the MSTs for different EN 
strengths. Fig.|4]upper panels show that as the EN mag- 
nitude increases, the steepness of the curve as a function 
of 0^0 is reduced for both tiow^m and Thi^iow Such 
changes in the MSTs serve to make the less favorable 
state more populated across a wide range of xq values. 




FIG. 4: (Color online) (a) tiow^m and (b) Thi^iow for various 
EN strengths and Tc = 10^. The population fraction Piow (c) 
and Phi (d) in the low and high states in steady state. 



To illustrate this effect, the lower panels of Fig. [4] show 
the probability of the system being in the low or hi state. 



calculated as Plow = ThiUiow 



/o 



hi-^low 



^low^hi) and 



Phi = 1 — Plow, respectively. Here one can see that as 
the EN magnitude is increased, not only does the abso- 
lute range of bistability expand but so does the range at 
which the population is macroscopically heterogeneous 
{e.g. 1 part in 100). The tails of these probabilities de- 
crease much more slowly than a system with only IN. 
Conclusions. Considered in the context of a population 
of cells, our analysis of a simple SRG model shows that 
EN is one of the primary drivers of phenotype switch- 
ing. Switching times can be lowered by multiple orders 
of magnitude and the mechanism of switching may not 
be strictly IN-driven, as previously assumed. If we inter- 
pret Xq in our model as an environmental input (e.^., the 
concentration of an inducer or antibiotic), then the pa- 
rameter range at which a cellular population will exhibit 
macroscopic levels of heterogeneity is greatly expanded 
by EN. Also, by showing how EN can modify the tails 
of bistable PDFs, our theory provides an interpretation 
for experimental observations of cells persisting in lowly 
populated phenotypes across unexpected conditions. 

Although this study was based on a simple SRG 
model, the general results apply to more complex ge- 
netic switches where EN is present in many kinetic rates. 
Stochastic models of cellular decision making will need to 
account for EN if they are to correctly recover switching 
times and trajectories. However, the major roadblock is 
the lack of experimental data regarding the properties 
of EN. It may be possible to use our theory to deconvo- 
lute the effects of IN and EN on switching from switching 
trajectories of individual cells subject to external fluctua- 
tions. We plan to explore such possibilities in the future. 
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